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We have developed an efficient order- iV real space Kubo approach for the calculation of the phonon 
conductivity which outperforms state-of-the-art alternative implementations based on the Green's 
function formalism. The method treats efficiently the time-dependent propagation of phonon wave 
packets in real space, and this dynamics is related to the calculation of the thermal conductance. 
Without loss of generality, we validate the accuracy of the method by comparing the calculated 
phonon mean free paths in disordered carbon nanotubes (isotope impurities) with other approaches, 
and further illustrate its upscalability by exploring the thermal conductance features in large width 
edge-disordered graphene nanoribbons (up to ~ 20nm), which is out of the reach of more con- 
ventional techniques. We show that edge-disorder is the most important scattering mechanism for 
phonons in graphene nanoribbons with realistic sizes and thermal conductance can be reduced by a 
factor of ~10. 
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I. INTRODUCTION 

In recent years, the understanding of phonon trans- 
port in carbon-based materials such as carbon nanotubes 
(CNTs) and graphene-based materials 2 has become par- 
ticularly important, both for fundamental studies of co- 
herent transport and also in view of novel applications. 
The thermal conductivity of suspended and supported 
single graphene layers has been found to be extremely 
high' 1 ' 4 owing to micrometer long phonon mean free 
paths (MFP). Such high thermal conductivity of two- 
dimensional graphene increases its potential for faster 
nano-electronic devices with less energy dissipation'. On 
the other hand, damaging a material like graphene could 
present interesting opportunities like achieving a high 
thermoelectric figure of merit and observation of An- 
derson localization of phonons. The question whether 
or not Anderson localization of acoustic phonons can 
be demonstrated unambiguously in disordered materials 
has been a long-standing problem in phonon physics, a 
phenomenon originating from the interference between 
multiple scattering paths was found ubiquitous in wave 
physics' 1 ''. Besides, disordered CNT-based bundles have 
been found to exhibit a tendency towards a thermal in- 
sulating regime'''''. Also, isotope disorder was shown to 
strongly impact on the high-energy phonon modes, re- 
sulting in very low mean free paths, in marked contrast 
with the genuine robustness of ballistic conduction for 
low-energy phonon modes 10 . As a result, the contribu- 
tion of quantum localization effects of high energy modes 
was found to be important but completely masked when 
considering the thermal conductance of the disordered 
material (being an integrated quantity over the entire 
phonon spectrum). Graphene nanoribbons (GNRs) of- 
fer an alternative to carbon nanotubes. By constructing 



heterostructures from pristine graphene mixed with dis- 
ordered (isotope impurities) GNRs 11 , or selective cover- 
age of hydrogen impurities 12 , it is possible to tune the 
resulting phonon transport capability. Another impor- 
tant source of disorder is provided by unavoidable ribbon 
edge irregularities, absent in CNTs which have already 
been shown to significantly reduce thermal conductance 
in small width GNRs 13 ' 14 . 

Exploring such possibilities, however, demand the de- 
velopment of efficient computational approaches which 
are able to tackle large scale (and realistic) simulations 
of material of interest. The Green's function (GF) meth- 
ods are able to incorporate microscopic details of the sys- 
tem, but they require matrix inversions which unavoid- 
ably limit the accessible size of the simulated systems. 
For the case of GNRs, although the system length can 
be upscaled without difficulty thanks to the decimation 
procedure, the computational cost becomes prohibitive 
for lateral sizes above 10 nm. 10,13 

In this Rapid Communication, using the real space 
Kubo (RSK) formalism we first demonstrate that a time- 
dependent phonon wave packet formalism can be con- 
nected to the calculation of the thermal conductance. 
After validating this numerical approach by comparing 
the obtained phonon MFPs in disordered CNTs (with 
isotope impurities) with previously computed ones by 
means of GF-based method 10 , we apply the new algo- 
rithm to large width GNRs, and focus on the impact of 
edge-disorder profiles. Scaling properties of phonon MFP 
and temperature-dependent thermal conductance are cal- 
culated as a function of edge-disorder strength and for 
lateral ribbon sizes accessible to today's state-of-the-art 
lithography. The broad generality of this method could 
offer a novel framework to explore other types of complex 
materials. 



freedom belongs. After some algebra, a becomes 



FIG. 1. (Color online) A short portion of an edge-disordered 
ZGNR with width w=17.04 nm, length L ~ 50 nm and dis- 
order density of 10%. In the calculation, the length of ZGNR 
is chosen to be 983.80nm, and periodic boundary condition is 
employed. 



II. COMPUTATIONAL PHONON TRANSPORT 
METHODOLOGY 

To investigate bulk quantum phonon transport in dis- 
ordered materials, the use of the Kubo formalism turns 
out to be the most natural and computationally efficient 
one. It has already been used for investigating ther- 
mal transport in disordered binary alloys or nanocrys- 
tralline silicon '. Similarly, many efforts have been 
made to explore time-dependent features of propagating 
phonon (or polaron) wave packets in complex quantum 
systems 1 '. Here a novel real space implementation of 
the Kubo formula for phonon propagation is given, es- 
tablishing a direct computational bridge between phonon 
dynamics and the thermal conductance. In comparison 
to the previous methods, we extract the dynamical in- 
formation from time evolution of the wave packet 16 and 
simulate Hamiltonian dynamics in place of Newtonian 
equations of motion 18 , so that only one initial condi- 
tion is required, i.e. the initial atomic displacements, 
and one does not need to calculate the atomic velocities 
in time. Additionally, by using the Lanczos technique 
which avoids any matrix inversion, a considerable com- 
putational efficiency gain is obtained, allowing the study 
of very large scale materials. The starting Hamiltonian 
which describes the phonon spectrum, taking only the 
harmonic interactions into account, is of the form 
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where Ui and pi are the displacement and momentum 
operators for the ith degree of freedom, Mi is the corre- 
sponding mass, and $ is the force constant tensor. Based 
on the linear response theory, the phonon conductivity a 
can be defined as VLT" 1 // d\ f™ dt(J x (-ih\)J x (t)) 16 . 

J x is the x component of the energy flux operator J, and 
it can be expressed as J x = l/2fi^\. (Aj — Xj)&ijUjVj, 
where ij is the velocity operator and A, is the equilib- 
rium position of the atom to which the ith degree of 



dw a W ^Tr{[A, D]S(oj 2 -D)[X, D]S(lo 2 -D)}, 
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(2) 
where fg is the Bose distribution function, and Z),j = 

&ij/ \jMi~Mj is the mass normalized dynamical matrix, 
Tr denotes the trace of the matrix, and S(u>2 — D) is the 
Dirac-delta operator. . The dynamical matrix can be 
obtained given the atomistic configuration of the system 
and it defines a Newtonian equation of motion which is 
of second order in time. Using the fact that the operators 
D and H have the same energy spectrum (Ti 2 D = H 2 ) 
and defining V x = [X,H]/ih, one can write the thermal 
conductance of a one-dimensional system as 
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(3) 
The thermal conductance can also be derived from the 
Landauer formalism or the nonequilibrium Green's 
function approach 20 as 

1 

K = 

2tt 



dwnw-^T(w), 



(4) 



with T{ui) being the phonon transmission function. 
Comparing the two formulas, the transmission function 
is defined as 
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Tr{V x S(huj - H)V x 5(huj - H)}, (5) 



which has the same form as the electron transmission 
function derived from the Kubo-Greenwood formula , 

Tei{E) = ^-Tr{V x S(E - H el )V x 5(E - H el )}. (6) 

This equivalence allows us to implement an order N al- 
gorithm related to the Lanczos method, which has been 
successfully applied to electron conduction in complex 
materials 21 . The Lanczos method has been previously 
employed in the past for calculating the vibrational den- 
sity of states of disordered systems 22 . The starting point 
of the RSK method is the transmission function as ex- 
pressed in Eq. (6). One can rewrite T(w) in terms of the 
diffusion coefficient, 
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where the diffusion coefficient T> has the form 

1 Tr{(A > (t)-A(0)) 2 <5(a; 2 - J D)} 



V(oj,t) 



Tr{S(u 2 -D)} 



(8) 



where X(t) is the position operator in the Heisenberg 
picture. In the diffusive regime D(ui,t) = T> xaasL (ui), so 
that the transmission function reduces to 

7» = -^Tr {5(uj 2 - D)}V max (u;), (9) 
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FIG. 3. (Color online) Transmission spectra for different 
lengths of ZGNR(80) with edge-disorder of 10% (solid) and 
15% (dashed). 



FIG. 2. (Color online) Elastic MFP for ZGNR of widths 
N z =20, 40, and 80 (4.26, 8.52, and 17.04 nm, respectively) 
with disorder density of 10%, and also for the N, =80 and 
15% disorder for comparison. Inset: Time dependent dif- 
fusion coefficients T>{ui,t) for three chosen frequencies. The 
same calculation time corresponds to different evolution time 
for different frequencies. 



where in the ballistic regime T>(u,t) = v 2 (uj)t, with 
v(ui) being the average group velocity over states 
with frequency u>. Since the number of channels 
is N ch (cj) w 2ujttTt {S{uj 2 -D)}v(u)/L, the phonon 
mean free path can be approximated by £{ui) = 
'Pmax (<•<->) /i>(w). By computing T>(u>,t), one can thus de- 
duce P max (w) and v(u)) and £(u>). Note that the numer- 
ator in Eq. (8) can be rewritten as Tr{[X,U(t)]^S(u} 2 — 
D)[X, U(t)]}, and approximated by Tr{[A, U{t)]H{uj 2 - 
D)[X,U(t)]} to the first order in perturbation the- 
ory, with U(t) = c~ iDt , and r = t/2u. The trace 
can be efficiently calculated through an average over 
a few random phase states of atomic displacements as 



N(xP\\X,U{t)V5{uj 2 - D)[X,U(t)]\i/}), N being the num- 
ber of degrees of freedom and the bra-ket corresponds to 
the local density of states (LDOS) associated with the 
vector [X,U(t)]\iP) which is calculated by using Cheby- 
shev expansion of U(t). LDOS is obtained by using 
the Lanczos method, and Tr {S(u 2 — D)} is also calcu- 
lated by averaging the LDOS of \ip) through the Lanczos 
method. 



III. RESULTS 

As a test case, we first consider a CNT(7,0) with 10.7% 
14 C impurities, which was studied in Ref. 10 using the 
GF method. From the saturation values of the time- 
dependent diffusion coefficients, we obtain the phonon 
MFP using the RSK scheme which compares very well 
with those obtained from the GF method 2 ' 5 . 

Next, we consider zigzag graphene nanoribbons 
(ZGNR) of different widths with edge-disorder. We use 



the fourth nearest neighbor force constants for building 
the dynamical matrices 24 . The ribbon widths are defined 
with the number of zigzag chains N z = 20, 40, and 80, 
and the relative amount of edge defects (removed carbon 
atoms at the edges) is chosen to be 10%, and additionally 
15% for N z = 80 (Fig. 1). The inset of Fig. 2 displays 
the evolution of the wave packet dynamics for different 
frequencies for ZGNR(80) with 10% edge-disorder.. The 
linear increase in T>{uj, t) at t > indicates ballistic trans- 
port at relatively short distances, whereas the decrease 
in T>(uj,t) is a signature of localization for this partic- 
ular frequency. The fact that T> decreases slowly sug- 
gests the localization length is large. The saturation of 
T>(uj,t) to a maximum value characterizes diffusive trans- 
port (c.f. Eq. 9). In Fig. 2, the decay of £(u>) with de- 
creasing the ribbon width is shown at a fixed disorder 
strength. This behavior can be rationalized with the fact 
that the scattering rate decreases with increasing width, 
a behavior previously derived for electron transport in 
both disordered CNTs and GNRs J '. One notes that, for 
low frequency modes, the MFP are several hundreds of 
nanometers, and due to large values of N c h the possibility 
to observe any onset of Anderson localization is jeopar- 
dized in the thermal conductance, as previously discussed 
for small diameter disordered carbon nanotubes 1 ". We 
then focus on ZGNR(80), and obtain the transmission 
according to T(u>) = A c h/(1 + L/£(tu)), assuming a dif- 
fusive regime, disregarding quantum interference effects 
and assuming minimum contact resistances. The result- 
ing frequency-dependent transmission function 7~(w) for 
different ribbon lengths is plotted in Fig. 3, while the 
thermal conductance k using Eq. (4) is shown in Fig. 4. 
The downscaling of T{u>) directly impacts on n which 
is found to be reduced by one order of magnitude for 
2 /um (at room temperature) compared with the ballistic 
case. Note that the calculation time directly determines 
the largest r that we can reach, the same r corresponds 
to different evolution time for different frequencies, since 
t = t/2cu. It takes longer calculation time for low fre- 
quency phonons to reach T> mllK . Here for u> < 70 cm -1 , 
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FIG. 4. (Color online) Temperature and length-dependent be- 
haviors of the ZGNR(80) thermal conductance (10% disorder 
(solid lines), 15% (dashed lines). 



time. To compute the contribution of these modes to k, 
a linear extrapolation for the transmission is used and 
the results are compared with those obtained from the 
GF method. Our analysis shows that this approximation 
causes an error of less than 1.5% for the thermal conduc- 
tance at room temperature. The conductivity of edge- 
disordered GNRs is obtained using a = kL/A, A being 
the cross section area of the ribbons which is taken as the 
interplane distance of graphite layers. At room tempera- 
ture, a = 1004 Wm-'K" 1 and a = 757 Wm" 1 ^ 1 with 
edge-disorder 10% and 15%, respectively for L =2 /xm. 
Comparing these values with the experimental values for 
suspended' (3000-5000 Win" 1 !^ 1 ) and supported l (600 
Wm -1 K -1 ) graphene, we conclude that edge-disorder is 
a very important source of scatterings not only in ultra- 
narrow ribbons, but also in GNRs as wide as at least 20 
nm. A comparison of our results with that of Ref. 13 
suggests that the effect of edge-disorder is reduced when 
increasing the width from 2 nm to 20 nm although it still 
plays an important role for large lateral sizes. The ratio 
of edge atoms affected from edge reconstructions to the 
total number of atoms decay inversely with the width of 



the GNR 11 , therefore edge profile disorder predominates 
over reconstruction effects with increased ribbon width. 



IV. CONCLUSION 

Conclusion.-A new real space and order N approach 
to compute the phonon wave packet propagation and 
the thermal conductance has been reported. Diffusion 
coefficients and MFPs can be extracted directly from 
the wave packet propagation. Its computational accu- 
racy and efficiency were demonstrated on disordered car- 
bon nanotubes and large width graphene nanoribbons 
respectively. A strong impact of smooth edge-disorder 
on the thermal conductance was found. Unlike edge 
reconstruction 1 , edge-disorder strongly suppresses ther- 
mal conduction not only for ultra- narrow GNRs, but 
also for realistically large ribbons. Phonons in edge- 
disordered GNRs, being scattering so effectively, pin- 
points towards good thermoelectric properties of large 
width GNRs. The applicability of this new method goes 
far beyond quasi-one dimensional systems and the area 
of carbon-based materials (nanotubes, graphene, etc.) 
studied here, and could be applied without difficulty to 
a wide range of other materials, including Boron-nitride- 
based materials 2 '' or silicon-based materials (nanowires, 
super lattices, etc.) 2 '. 
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